%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DTC_CD_Normalized_Baseline_fun.m; used by DTC_Normalized_Baseline_scr.m
% Edited by David Stern 17 March 2019 for final submission to JAERE.
% The baseline version finds the two normalization constants.
% This counterfactual simulation uses the numerical values from that program
% as its normalization constants.

%Export variable results to and import parameter values from script file
function [yratio, Y, Ym,Ys,Ypc, pratio,Nmt, Nst, Bm, Bs, B, n, Em, Es, es, em, eratio, EI, SshrY, SshrL, SshrE, SshrEexp, gYpc,gNm,gNs,Lratio,L] = DTC_CD_Normalized_Counterfactual_fun(gamma,Nm0,etam,etas,sigma,beta,alpha,Ns0,Embar,esbar,nu,T,p0,theta,L0,Ysb,Ymb)

h=(1-nu)/nu;
GAMMA=gamma/(1-gamma);


Lm=zeros(T,1);
Ls=zeros(T,1);
wage=zeros(T,1);
es=zeros(T,1);
em=zeros(T,1);
ps=zeros(T,1);
pm=zeros(T,1);
Profits=zeros(T,1);
Profitm=zeros(T,1);
Es=zeros(T,1);
Em=zeros(T,1);
Ys=zeros(T,1);
Ym=zeros(T,1);
Y=zeros(T,1);
Ypc=zeros(T,1);
Nmt=zeros(T+1,1); %NOTE t+1 for Nm means t for most other variables!
Nst=zeros(T+1,1); %NOTE t+1 for Nm means t for most other variables!
Nmt(1)=Nm0;
Nst(1)=Ns0;
Bm=zeros(T,1);
Bs=zeros(T,1);
L=zeros(T,1);
gNm=zeros(T,1);
gNs=zeros(T,1);
Lratio=zeros(T,1);
pratio=zeros(T+1,1);
pratio(1,1)=p0;

for t=1:T
    %population growth Equation (47)
    %L(t,1)=L0*(9.7/(1+exp(-log(81)/267*(1540+(t)*20-1530)))+47.4/(1+exp(-log(81)/171*(1540+(t)*20-1870))))/(9.7/(1+exp(-log(81)/267*(1560-1530)))+47.4/(1+exp(-log(81)/171*(1560-1870))));
    L(t)=1;
    %L(t)=exp(((t)-1)*.0001);
    %Fixed quantity of wood and price of coal.
    Em(t,1) = Embar; %Em(t) needed not for calc, but to report E, Emexp, etc
    es(t,1) = esbar; %es(t) needed not for calc, but to report Esexp, etc
    
    %define functions in Nm, Ns and p
    ps_fun=@(p) (gamma^sigma*p^(1-sigma)+(1-gamma)^sigma)^(1/(sigma-1)); %A1
    pm_fun=@(p) p*ps_fun(p); %A2
    Ls_fun=@(p) L(t,1)/(1+((gamma/(1-gamma))^sigma)*p^(1-sigma)); %A3
    Lm_fun=@(p) L(t,1)-Ls_fun(p); %A4
    Es_fun=@(Ns,p) (beta*es(t,1)/(alpha*(Ns+theta*Nst(t,1))))^((beta-1)/(1-alpha-beta))*Ls_fun(p)*(ps_fun(p)/Ysb)^(1/(1-alpha-beta)); %19 (in Appendix A1)
    Ym_fun=@(Nm,p) (1/beta)*(Nm+theta*Nmt(t,1))*(pm_fun(p)/Ymb)^(beta/(1-beta))*Em(t,1)^(alpha/(1-beta))*Lm_fun(p)^((1-alpha-beta)/(1-beta)); %17 (in Appendix A1)
    Ys_fun=@(Ns,p) (1/beta)*(Ns+theta*Nst(t,1))*(ps_fun(p)/Ysb)^(beta/(1-beta))*Es_fun(Ns,p)^(alpha/(1-beta))*Ls_fun(p)^((1-alpha-beta)/(1-beta)); %18 (in Appendix A1)
    Profitm_fun=@(Nm,p) (1-beta)*(pm_fun(p)/Ymb)^(1/(1-beta))*Lm_fun(p)^((1-alpha-beta)/(1-beta))*Em(t,1)^(alpha/(1-beta)); %14
    Profits_fun=@(Ns,p) (1-beta)*(ps_fun(p)/Ysb)^(1/(1-beta))*Ls_fun(p)^((1-alpha-beta)/(1-beta))*Es_fun(Ns,p)^(alpha/(1-beta)); %14


    %equation system for Nm, Ns and p: (A11)=(19), (A12)=(20), (A7)=(18)
     DTC_fun=@(Nm,Ns,p) [((Nm-Nmt(t,1))^h)/(nu*(etam^(1/nu)))-Profitm_fun(Nm,p);... %29
                           ((Ns-Nst(t,1))^h)/(nu*(etas^(1/nu)))-Profits_fun(Ns,p);... %30
                           p-(gamma/(1-gamma))*((Ym_fun(Nm,p)/Ymb)/(Ys_fun(Ns,p)/Ysb))^(-1/sigma)]; %28
         
     % solver
     options=optimset('TolX',1e-8,'TolFun',1e-8);
     DTC=fsolve(@(DTC) DTC_fun(DTC(1),DTC(2),DTC(3)),[Nmt(t,1),Nst(t,1),pratio(t,1)],options); 
     % Transfers results to variables
     Nmt(t+1,1)=DTC(1);
     Nst(t+1,1)=DTC(2);
     pratio(t+1,1)=DTC(3);
 
    
     % Formulae for display/output:
     ps(t,1)=ps_fun(pratio(t+1,1));
     pm(t,1)=pm_fun(pratio(t+1,1));
     Lm(t,1)=Lm_fun(pratio(t+1,1));
     Ls(t,1)=Ls_fun(pratio(t+1,1));
     Lratio(t,1)=Lm(t,1)/L(t,1);
     Es(t,1)=Es_fun(Nst(t+1),pratio(t+1,1));
    
     Profitm(t,1)=Profitm_fun(Nmt(t+1,1),pratio(t+1,1));
     Profits(t,1)=Profits_fun(Nst(t+1,1),pratio(t+1,1));
     Ym(t,1)=Ym_fun(Nmt(t+1,1),pratio(t+1,1));
     Ys(t,1)=Ys_fun(Nst(t+1,1),pratio(t+1,1));
     Y(t,1)=(gamma*((Ym(t,1)/Ymb)^((sigma-1)/sigma))+(1-gamma)*((Ys(t,1)/Ysb)^((sigma-1)/sigma)))^(sigma/(sigma-1)); % (1)
     Ypc(t,1)=Y(t,1)./L(t,1);
     em(t,1)=alpha*pm(t,1)*Ym(t,1)/(Em(t,1)*Ymb); % (7)
     gNm(t,1) = (Nmt(t+1,1)-Nmt(t,1))/Nmt(t,1);
     gNs(t,1) = (Nst(t+1,1)-Nst(t,1))/Nst(t,1);
     wage(t,1)=(1-alpha-beta)*pm(t,1)*Ym(t,1)/(Ymb*Lm(t,1)); % (10)
     Bm(t,1)=Nmt(t+1,1)+theta*Nmt(t,1);
     Bs(t,1)=Nst(t+1,1)+theta*Nst(t,1);

end

% Shares and expenditure data
Svalue=ps.*Ys/Ysb;
SshrY=Svalue./Y; %same as Sshare
SshrL=Ls./L;
E=Em+Es;
Emexp=em.*Em;
Esexp=es.*Es;
Eexp=Emexp+Esexp;
SshrE=Es./E;
SshrEexp=Esexp./Eexp;

%Output ratio
yratio=Ym./Ys;
%Energy intensity
EI=(Em+Es)./Y;
%Annual Economic Growth Rate

gYpc=(((diff(Ypc)./Ypc(1:end-1)+1)).^0.05-1)*100;
%Energy price ratio
eratio = em./es;
%Backwardness
B=Bm./Bs;
%Direction of technical change
n=GAMMA^(1/h)*B.^(-(1+h)/h).*yratio.^((sigma-1)/(h*sigma));



